路遙知碼力，日久練成精- 只要在程式之路鑽研的夠深，便能夠充分發揮程式碼的力量; 練習的日子夠久，便能夠練成寫出精簡代碼的能力。

經歷了十五天python精華語法介紹的前奏後，
終於要來做第一個程式小專題囉。

古希臘人相信每個數字都有自己代表的意義，
例如1表示創造，2表示女性，3代表男性，
2+3變成的5代表結婚，
接下來則是特別的6，
6的因數有1, 2, 3, 6，
6除了自己之外的因數和恰好是1+2+3=6，
因此6被視為是一個「完美數」。

今天的主題主要是探討兩種特別的數對:「友好數」、「婚約數」。
在以前計算機、電腦還不發達的年代，
要找到諸如此類特別的數字，
只能靠數學家用大量的數學技巧去慢慢發現。

今日，藉著電腦高速的運算速度，
即便只有中學的數學知識，
你我都有能力找出這些特別的數，
以滿足人類自古以來的好奇心

完美數在數學界的名聲比較大，
數學上已經有大量研究了，
今天我們主要來找友好數(amicable pair)、婚約數（betrothed numbers）。

友好數、婚約數

首先給個定義:
友好數就是兩個正整數，
彼此除了自身之外的正因數之和與對方相等。
例如 (220, 284) 是一對最小的友好數。
220 自身外的正因數和為1+2+4+5+10+11+20+22+44+55+110=284
284 自身外的正因數和為1+2+4+71+142=220
所以友好數就像一對好朋友一樣，
彷彿一種「你中有我，我中有你」的感情。

婚約數的定義跟友好數非常相似，
婚約數是彼此除了1與自身之外的正因數之和與對方相等。
例如 (48, 75) 是一對最小的婚約數。
48除了1和本身的正因數相加是：2+3+4+6+8+12+16+24=75，
75除了1和本身的正因數相加是：3+5+15+25=48。

第一步: 求出整數的正因數和

希望大家還記得什麼是因數，
n的因數就是可以整除n的數字。
我們運用之前學過的列表生成式，
達到方便求和的效果。
以下程式可以求出除了自身外的正因數和。

def divisorsum(n):
    return sum([i for i in range(1,n) if n%i==0])

相信這是最直覺的方法，
我們從1開始搜索，
利用 n%i==0 這個判斷式判斷是否為因數，
如果是因數就添加到列表中，
最後用sum()函數求和。

第二步: 雖然非主軸但順便教大家求完美數

有了求正因數和的函數，
我們可以很容易判斷一個數是否為完美數。

def isPerfect(n):
    return divisorsum(n) == n

我們定義一個函數isPerfect(n)，
判斷一個數是不是完美數，
即判斷一個數自身外的正因數和是否等於自己。
加入以下程式找出小於10000的完美數

for i in range(1,10000):
    if isPerfect(i):
        print(i)

結果:
6
28
496
8128
你可以參考維基百科，
確認我們找到的完全數正確無誤。

第三步: 更有效率的求出整數的正因數和

求正因數和的程式碼很精簡，
但目前我們仍然不太滿意，
當我們將上述程式改為尋找小於十萬的完美數:

for i in range(1,100000):
    if isPerfect(i):
        print(i)

你可以試著在你的電腦執行看看，
執行時間很可能會跑個幾十秒。
因為我們的計算因數和的函數是沒效率的。
我們如何找到n的所有因數呢?
原本我們讓程式去檢查所有1~n-1的數字是否可以整除它
但實際上，我們只要搜索到根號n就可以找到所有因數。

對n來說，
我們把小於根號n 的數字稱為「小數字」，
大於根號n的數字稱為「大數字」。
例如根號30大約是5點多，

30 = 1*30
   = 2*15
   = 3*10
   = 5*6

每個分解都必然是一個「小數字」乘上「大數字」。
你硬是讓兩個「大數字」相乘，
乘出來的數字就太大了。
(兩個大於根號n的數相乘一定大於n)

所以說，我們只要能夠找到小於根號30的因數「1, 2, 3, 5」，
便能透過除法去算出大於根號30的因數「30, 15, 10, 6」。

先教大家如何在python開根號，
開根號在數學上相當於「0.5次方」，
例如:

>>> 2 ** 0.5
1.4142135623730951
>>> 4 ** 0.5
2.0
>>> 9 ** 0.5
3.0

但是注意用次方運算出來的數是帶有小數點的浮點數，
我們可以取int()把小數點去掉，例如:

>>> 30 ** 0.5
5.477225575051661
>>> int(30 ** 0.5)
5

根據上述想法，
我們可以把求正因數和的計算式改寫如下:

def divisorsum(n):
    return sum([sum([i,n//i]) for i in range(1,int(n**0.5)+1) if n%i==0])-n

看這段代碼可能不容易理解，
我們試著印出列表中的內容看看:
(假設n是30)

print([[i,30//i] for i in range(1,int(30**0.5)+1) if 30%i==0])

結果為 [[1, 30], [2, 15], [3, 10], [5, 6]]
這不就是上例分解30為兩個整數乘積的狀況嗎?
我們把列表中的數字全部相加就是因數和了。
但是注意若我們把30改為9 (9本身是完全平方數)，
結果為 [[1, 9], [3, 3]]，
「根號9」這個因數會被重複計算，
簡潔起見，我們把列表內的列表改為集合以去除重複。

print([{i,9//i} for i in range(1,int(9**0.5)+1) if 9%i==0])

結果變為[{1, 9}, {3}]。

因此原程式改為:

def divisorsum(n):
    return sum([sum({i,n//i}) for i in range(1,int(n**0.5)+1) if n%i==0])-n

注意因為我們這邊要計算除了「自身」外的正因數和，
計算的結果還要再減去n。

第四步: 求出小於十萬的友好數數對

先做個小測試，
我們用簡單的列表生成式算出1~6除了自身外的正因數之和:

print([(i,divisorsum(i)) for i in range(1,7)])

結果為:
[(1, 0), (2, 1), (3, 1), (4, 3), (5, 1), (6, 6)]
意思是
1(除了自身)的正因數和=0，
2(除了自身)的正因數和=1，
…
6(除了自身)的正因數和=6。
大家可以自行驗算。

還記得什麼是友好數嗎?
友好數的條件是自身外的正因數和等於彼此，
例如(220, 284) 是一對最小的友好數，
divisorsum(220)=284，divisorsum(284)=220。
也就是說，若n滿足divisorsum(divisorsum(n))==n的條件，
則(n, divisorsum(n))便是一組友好數。

噢，對了，divisorsum這個函數的名字有點長，
這邊故且以DS這個名字代替。
改寫如下:

def DS(n):
    return sum([sum({i,n//i}) for i in range(1,int(n**0.5)+1) if n%i==0])-n

另外我們定義amicablePair(k)函數，
計算小於k的友好數，
判斷條件就用 DS(DS(i))==i 來判斷，如下:

def amicablePair(k):
    return [(i,DS(i)) for i in range(1,k+1) if DS(DS(i))==i]

我們先試試看效果如何，
測試找300以下的友好數:

print(amicablePair(300))

結果為[(6, 6), (28, 28), (220, 284), (284, 220)]

檢查答案有誤，
(6, 6)因為是同一個數，
並不是我們想要的答案。
另外(220, 284), (284, 220)其實算同一組數對，
被重複計算了。
我們還要再加上一條判斷條件i< DS(i)才完整，
調整如下:

def amicablePair(k):
    return [(i,DS(i)) for i in range(1,k+1) if i< DS(i) and DS(DS(i))==i]

重新測試結果:

print(amicablePair(300))

結果為[(220, 284)]，
這次已經能正確找到一組數對了。
試試找出小於十萬的友好數吧:

print(amicablePair(100000))

結果為[(220, 284), (1184, 1210), (2620, 2924), (5020, 5564), (6232, 6368), (10744, 10856), (12285, 14595), (17296, 18416), (63020, 76084), (66928, 66992), (67095, 71145), (69615, 87633), (79750, 88730)]。

第五步: 舉一反三- 求出小於十萬的婚約數數對

婚約數與友好數定義很類似，
婚約數是彼此除了1與自身之外的正因數之和與對方相等。
我們定義函數DS2(n) 計算n除了1與自身之外的正因數之和:

def DS2(n):
    return sum([sum({i,n//i}) for i in range(1,int(n**0.5)+1) if n%i==0])-n-1

只要剛剛的講解你都有跟上的話，
計算婚約數的函數也只要微調一下即可:

def betrothedPair(k):
    return [(i,DS2(i)) for i in range(1,k+1) if i< DS2(i) and DS2(DS2(i))==i]

再附上小於十萬的婚約數結果:

print(betrothedPair(100000))

結果: [(48, 75), (140, 195), (1050, 1925), (1575, 1648), (2024, 2295), (5775, 6128), (8892, 16587), (9504, 20735), (62744, 75495)]

第六步: 大功告成，總結所有程式碼

我們將今天講解的所有程式碼統整在一起，
並順便計算出小於十萬的完美數。
完整的程式碼如下:

def DS(n):
    return sum([sum({i,n//i}) for i in range(1,int(n**0.5)+1) if n%i==0])-n

def DS2(n):
    return sum([sum({i,n//i}) for i in range(1,int(n**0.5)+1) if n%i==0])-n-1

def findPerfect(k):
    return [i for i in range(1,k+1) if DS(i) == i]

def amicablePair(k): #求小於k的友好數對
    return [(i,DS(i)) for i in range(1,k+1) if i< DS(i) and DS(DS(i))==i]

def betrothedPair(k): #求小於k的婚約數對
    return [(i,DS2(i)) for i in range(1,k+1) if i< DS2(i) and DS2(DS2(i))==i]
    
print(findPerfect(100000))
print(amicablePair(100000))
print(betrothedPair(100000))

結果:
[6, 28, 496, 8128]
[(220, 284), (1184, 1210), (2620, 2924), (5020, 5564), (6232, 6368), (10744, 10856), (12285, 14595), (17296, 18416), (63020, 76084), (66928, 66992), (67095, 71145), (69615, 87633), (79750, 88730)]
[(48, 75), (140, 195), (1050, 1925), (1575, 1648), (2024, 2295), (5775, 6128), (8892, 16587), (9504, 20735), (62744, 75495)]

名字彩蛋

結果解完了，不知讀者們有發現友好數及婚約數的命名由來嗎?
目前已知的友好數對全部都是「兩個奇數」或同為「兩個偶數」的組合，
而目前已知的婚約數則全部都是「一奇一偶」的組合。
或許因為這樣的特性(相同奇偶性 => 同性; 相異奇偶性 => 異性)，
使得這樣的數對各自稱為友好數及婚約數呢。